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Abstract 

In this paper we perform an analytical and numerical study of Extreme Value distri- 
butions in discrete dynamical systems. In this setting, recent works have shown how 
to get a statistics of extremes in agreement with the classical Extreme Value Theory. 
We pursue these investigations by giving analytical expressions of Extreme Value dis- 
tribution parameters for maps that have an absolutely continuous invariant measure. 
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We compare these analytical results with numerical experiments in which we study 
the convergence to limiting distributions using the so called block-maxima approach, 
pointing out in which cases we obtain robust estimation of parameters. In regular 
maps for which mixing properties do not hold, we show that the fitting procedure 
to the classical Extreme Value Distribution fails, as expected. However, we obtain 
an empirical distribution that can be explained starting from a different observable 



function for which Nicolis et al. 120061 have found analytical results. 



1 Introduction 



Extreme Value Theory (EVT) was first developed by Fisher and Tippett |1928 



and formalized by |Gnedenko |1943 which showed that the distribution of the 



block-maxima of a sample of independent identically distributed (i.i.d) vari- 
ables converges to a member of the so-called Extreme Value (EV) distribution. 
It arises from the study of stochastical series that is of great interest in differ- 
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1989 . All these events have 


a relevant impact 



on socioeconomic activities and it is crucial to find a way to understand and, if 



possible, forecast them Hallerberg and Kantz 2008 , Kantz et al. 2006 



The attention of the scientific community to the problem of modeling extreme 
values is growing. An extensive account of recent results and relevant applica- 

Such an interest is mainly due to the fact 



Ghil et al. 2011 



tions is given in 

that this theory is also important in defining risk factor in a wide class of appli- 
cations such as the modeling of financial risk after the significant instabilities in 
financial markets worldwide 



et al. 1999 



Gilli andKellezi 2006 , |Longin 2000 



tins and Stedinger 



Embrechts 

the analysis of seismic and hydrological risk Burton 1979 , |Mar- 
2000 1 . Even if the probability of extreme events decreases 



with their magnitude, the damage that they may bring increases rapidly with 



the magnitude as does the cost of protection against them Nicolis et al. 2006 1. 
From a theoretical point of view, extreme values represent extreme fluctuations 
of a system. Very recently, many authors have shown clearly how the statis- 
tics of global observables in correlated systems can be related to EV statistics 
jDahlstedt and Jensen 2001 



Bertin[ |2005| . |Clusel and Bertinj (2008 have 



shown how to connect fluctuations of global additive quantities, like total en- 
ergy or magnetization , by statistics of sums of random variables in such a way 
that it is possible to identify a class of random variables whose sum follows an 
extreme value distributions. 

The so called block-maxima approach is widely used in EVT since it represents 
a very natural way to look at extremes. It consists of dividing the data series 
of some observable into bins of equal length and selecting the maximum (or 
the minimum) value in each of them Coles et al. 1999 . When dealing with 
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climatological or financial data, since we usually have limited data-set, the main 
problem in applying EVT is related to the choice of a sufficiently large statistics 
of extremes provided that each bin contains a suitable number of observations. 
Therefore a smart balance between number of maxima and observations per 
bin is needed 



et al. 2005 



Felici et al.[|2007a| , |Katz and Brown| |1992] , |Katz| |1999| , [Katz 



Recently a number of alternative approaches have been studied. One consists in 
looking at exceedance over high thresholds rather than maxima over fixed time 
periods. While the idea of looking at extreme value problems from this point of 
view is very old, the development of a modern theory has started with |Todorovic| 



and Zelenhasic 1970 that have proposed the so called Peaks Over Threshold 
approach. At the same time there was a mathematical development of proce- 
dures based on a certain number of extreme order statistics IPickands IIII 11975 



Hill 1975 and the Generalized Pareto distribution for excesses over thresholds 



Smith 1984 , Davison 1984 , Davison and Smith 1990 



Since dynamical systems theory can be used to understand features of phys- 
ical systems like climate and forecast financial behaviors, many authors have 
studied how to extend EVT to these field. When dealing with dynamical sys- 
tems we have to know what kind of properties (i.e. stability, degree of mixing, 
correlations decay) are related to Gnedenko's hypotheses and also which observ- 
ables we must consider in order to obtain an EV distribution. Furthermore, even 
if the convergence is achieved, we should evaluate how fast it is depending on 
all parameters and properties used. Empirical studies show that in some cases 
a dynamical observable obeys to the extreme value statistics even if the con- 



2007 , Vitolo et al. 


2009b , Vitolo et al. 


2009a. J 


et al. 


1995] and more recently 


Nicolis et al. 


2006 



For example, Balakrishnan 



and Haiman 2003 have 



shown that for regular orbits of dynamical systems we don't expect to find con- 
vergence to EV distribution. 

The first rigorous mathematical approach to extreme value theory in dynamical 
systems goes back to the pioneer paper by P. Collet in 2001 |Collet 2001 . Col- 
let got the Gumbel Extreme Value Law (see below) for certain one-dimensional 
non-uniformly hyperbolic maps which admit an absolutely continuous invariant 
measure and exhibit exponential decay of correlations. Collet's approach used 



Young towers Young 1999 , |Young 1998 and his suggestion was successively 



applied to other systems. Before quoting them, we would like to point out that 
Collet was able to establish a few conditions (usually called D and D') and 
which have been introduced by Leadbetter |Leadbetter et alT 



1983 with the 



aim to associate to the stationary stochastic process given by the dynamical 
system, a new stationary independent sequence which enjoyed one of the classi- 
cal three extreme value laws, and this law could be pulled back to the original 
dynamical sequence. Conditions D and D' require a sort of independence of 
the stochastic dynamical sequence in terms of uniform mixing condition on the 
distribution functions. Condition D was successively improved by Freitas and 
Freitas iFreitas and Freitas 2008a , in the sense that they introduced a new 
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condition, called £>2, which is weaker than D and that could be checked directly 
by estimating the rate of decay of correlations for Holder observables Q We 
notice that conditions D2 and D' allow immediately to get Extreme Value Laws 
for absolutely continuous invariant measures for uniformly one-dimensional ex- 
panding dynamical systems: this is the case for instance of the 1-D maps with 
constant density studied in Sect. 3 below. Another interesting issue of Collet's 
paper was the choice of the observables g's whose values along the orbit of the 
dynamical systems constitute the sequence of events upon which we successively 
search for the partial maximum. Collet considered a function <?(dist(x, £)) of the 
distance with respect to a given point £, with the aim that g achieves a global 
maximum at almost all points £ in the phase space; for example g{x) = — log a;. 
Using a different g, Freitas and Freitas Moreira Freitas and Freitas 2008 were 



able to get the Weibull law for the family of quadratic maps with the Benedicks- 
Carlesson parameters and for ( taken as the critical point or the critical value, 
so improving the previous results by Collet who did not keep such values in his 
set of full measure. 



The latter paper Moreira Freitas and Freitas 2008 strongly relies on condition 
D2; this condition has also been invoked to establish the extreme value laws 
on towers which model dynamical systems with stable foliations (hyperbolic 
billiards, Lozi maps, Henon diffeomorphisms, Lorenz maps and flows). This 
is the content of the paper by Gupta, Holland and Nicol |Gupta et al. 2009] . 
We point out that the observable g was taken in one of three different classes 
91)92,53, see Sect. 2 below, each one being again a function of the distance 
with respect to a given point £. The choice of these particular forms for the 
g's is just to fit with the necessary and sufficient condition on the tail of the 
distribution function F(u), see next section, in order to exist a non-degenerate 



limit distribution for the partial maxima Freitas et al. 2009 , Holland et al 



2008 . The paper Gupta et al. 2009 also covers the easier case of uniformly 
hyperbolic diffeomorphisms, for instance the Arnold Cat map which we studied 
in Sect. 3.2. 

Another major step in this field was achieved by establishing a connection be- 
tween the extreme value laws and the statistics of first return and hitting times, 



see the papers by Freitas, Freitas and Todd |Freitas et al. 2009 , |Freitas et al. 



2010b . They showed in particular that for dynamical systems preserving an 



absolutely continuous invariant measure or a singular continuous invariant mea- 



1 We briefly state here the two conditions, we defer to the next section for more details 
about the quantities introduced. If X n ,n > is a stochastic process, we define Mj i = 
{Xj , Xj+i , ■ ■ ■ , Xjj^i } and we put Mo, m = M m ■ Moreover we set a m and b m two normalising 
sequences and u m = x/a m + b m , where x is a real number, cf. next section for the meaning 
of these variables. The condition D2(u m ) holds for the sequence J(.m if for any integer l,t,m 
we have |^(Xq > u m ,M t ; < u m ) — u(Xq > u m )i>(M t i < u m )\ < 7(m,t), where r y(m,t) is 
non-increasing in t for each m and m7(m,t m ) — > as m — ¥ 00 for some sequence t m = o(m), 
t m — > 00. 

We say condition D'(u m ) holds for the sequence X m if lim^^oo limsup m m^J™!*' u (Xq > 
u m ,Xj > u m ) = 0. Whenever the process is given by the iteration of a dynamical systems, 
the previous two conditions could also be formulated in terms of decay of correlation integrals, 
see |Preitas and Freitas] [2008b| , |Gupta| |2010| . 
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sure v , the existence of an exponential hitting time statistics on balls around 
v almost any point £ implies the existence of extreme value laws for one of the 
observables of type g^i — 1,2,3 described above. The converse is also true, 
namely if we have an extreme value law which applies to the observables of type 
gi,i = 1,2,3 achieving a maximum at £, then we have exponential hitting time 
statistics to balls with center £. Recently these results have been generalized 



to local returns around balls centered at periodic points Freitas et al. 2010a 



We would like to point out that the equivalence between extreme values laws 
and the hitting time statistics allowed to prove the former for broad classes of 
systems for which the statistics of recurrence were known, for instance for ex- 
panding maps in higher dimension. 



In this work we consider a few aspects of the extreme value theory applied 
to dynamical systems throughout both analytical results and numerical exper- 
iments. In particular we analyse the convergence to EV limiting distributions 
pointing out how robust arc parameters estimations. Furthermore, we check 
the consistency of block-maxima approach highlighting deviations from theo- 
retical expected behavior depending on the number of maxima and number of 
block-observation. To perform our analysis we use low dimensional maps with 
different properties: mixing maps in which we expect to find convergence to EV 
distributions and regular maps where the convergence is not ensured. 
The work is organised as follow: in section 2 we briefly recall methods and 
results of EVT for independent and identical distributed (i.i.d.) variables and 
dynamical systems. In section 3 we explicitly compute theoretical expected dis- 
tributions parameter in respect to the observable functions of type gi, i = 1,2,3 
for map that have constant density measure. Numerical experiments on low di- 
mensional maps are presented. In section 4 we show that it is possible to derive 
an asymptotic expression of normalising sequences when the density measure 
is not constant. As an example we derive the explicit expressions for the Lo- 
gistic map. Eventually, in section 5 we repeat the experiment for regular maps 
showing that extreme values laws do not follow from numerical experiments. 



2 Background on EVT 



Gnedenko] |1943| studied the convergence of maxima of i.i.d. variables 



Xq, X\,..X 

with cumulative distribution (cdf) F{x) of the form: 

F(x) = P{a m {M m - b m ) < x} 

Where a m and b m are normalising sequences and M rn = max{Vo, X±, V m _i}. 
It may be rewritten as F(u m ) = P{M m < u m } where u m — x/a m + b m . Such 
types of normalising sequences converge to one of the three type of Extreme 
Value (EV) distribution if necessary and sufficient conditions on parent distri- 
bution of Xi variables are satisfied Leadbetter et al. 1983 . EV distributions 
include the following three families: 
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• Gumbel distribution (type 1): 

F(x) = exp{— e _a: } x£ 

• Frechet distribution (type 2): 



F(x) = 
F(x) = exp {- 



.a/. 



x < 

£} x>o 



(i) 



(2) 



Weibull distribution (type 3): 



F(x) = exp(-(-x) x<0 
F(x) = l x>0 



(3) 



Let us define the right endpoint xf of a distribution function F{x) as: 

x F = sup{a; : F(x) < 1} (4) 

then, it is possible to compute normalising sequences a m and b m using the 
following corollary of Gnedenko's theorem : 

Corollary (Gnedenko): The normalizing sequences a rn and b m in the conver- 
gence of normalized maxima P{a m {M rn — b m ) < x} — > F(x) may be taken (in 
order of increasing complexity) as: 



• Type 1: 

• Type 2: 

• Type 3: 
where 



7m 1 , b m = Q or b m = c- to" 
(x F -J m )~ l , b m = x F ; 



= F _1 (l - 1/to) = M{x;F{x) > 1 - 1/m} 



G(t) = 



XF 1 - F(u) 



du, t < xf 



(5) 



(6) 



l-F(t) 

and c G K is a constant. It is important to remark that the choice of 
normalising sequences is not unique |Leadbetter et ah 1983 . For example for 
b m of type 2 distribution it is possible to choose either b m = or b m = c ■ mT^ . 
In particular, we will use the last one since it is a more general choice that ensure 
the convergence for a much broader class of initial distributions |Beirlant 2004 



Instead of Gnedenko's approach it is possible to fit unnormalized data directly 
to a single family of generalized distribution called GEV distribution with cdf: 



F G {x]fi,<7,t) = exp < - 



x — fl 



(7) 
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which holds for 1 + £(x — [i)ja > 0, using \i G E (location parameter) and 
cr > (scale parameter) as scaling constants in place of b m , and a m |Pickands III 
1968| . £ G K is the shape parameter also called the tail index: when £ — > 0, 
the distribution corresponds to a Gumbel type. When the index is negative, it 
corresponds to a Weibull; when the index is positive, it corresponds to a Frechet. 

In order to adapt the extreme value theory to dynamical systems, we will 
consider the stationary stochastic process Xo,Xi, ... given by: 

X m (x) = g(dist(f m (x), 0) Vm G N (8) 

where 'dist' is a Ricmannian metric on Q. £ is a given point and g is an 
observable function, and whose partial maximum is defined as: 

M m = max{X , . . . , X m _ J (9) 

The probability measure will be here an invariant measure v for the dy- 
namical system. As we anticipated in the Introduction, we will use three types 
of observables gi,i = 1,2,3, suitable to obtain one of the three types of EV 
distribution for normalised maxima: 

ffl (a;) = -log(diBt(x,0) (10) 
g 2 (x)=di S t(x,()- 1/a (11) 

g 3 (x) = C- dist (x,0 1/a (12) 

where C is a constant and a > G K. 
These three type of functions are representative of broader classes which are 
defined, for instance, throughout equations (1.11) to (1.13) in Freitas et al. 
[2009]; we now explain the reasons and the meaning of these choices. First of all 
these functions have in common the following properties: (i) they are defined 
on the positive semi-axis [0, oo] with values into R U {+oo}; (ii) is a global 
maximum, possibly equal to +oo; (iii) they are a strictly decreasing bijection in 
a neighborhood V of with image W. Then we consider three types of behavior 
which generalize the previous specific choices: 

Type 1: there is a strictly positive function p : W — > E such that Vy G E we 
have 

lim ^(s + ypjs)) =e _ y 
*->fli(o) g x (s) 
Type. 2: 32(0) = +00 and there exists /3 > such that Wy > we have 

S ^°° 92 ( s ) 
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Type 3: 173(0) = D < +00 and there exists 7 > such that Vy > we have 

^Og-l( D - S ) 



The Gnedenko corollary says that the different kinds of extreme value laws 
are determined by the distribution of F(u) = v(Xq < u) and by the right end- 
point of F, xp. We will see in the next section that the local invertibility of 
gi,i — 1,2,3 in the neighborhood of together with the Lebesgue's differenti- 
ation theorem (which basically says that whenever the measure v is absolutely 
continuous with respect to Lebesgue with density p, the the measure of a ball 
B$(xq) of radius d centered around almost any point xq scales like 8p(xo)), allow 
us to compute the tail of F, in fact we have 

l-F(u)~p(C)|B 9 -i( u )(0l, 

where g is any of the three types of functions introduced in (10) to (12) and \A\ 
denotes the diameter of the set A. As we said above the tail of F determines 
the three limit laws for partial maximum of i.i.d. sequences. In particular 



Th. 1.6.2. in Leadbetter et al. 1983 specifies what kind of conditions the 



distribution function F must verify to get one specific law: the above type 1,2,3 
assumptions are just the translation in terms of the shape of <?.; of the conditions 
on the tail of F. 



3 Distribution of Extremes in mixing maps with constant 
density measure 

Our goal is to use a block-maxima approach and fit our unnormalised data 
to a GEV distribution; for that it will be necessary to find a linkage among 
a m , b m , [i and a. At this regard we will use Gnedenko's corollary to compute 
normalising sequences showing that they correspond to the parameter we obtain 
fitting directly data to GEV distribution. 

We derive the correct expression for mixing maps with constant density 
measure and the asymptotic behavior for logistic map that is a case of non- 
constant density measure. 



3.1 Asymptotic sequences 

In this section we will consider the case of uniformly hyperbolic maps which 
preserve the Lebesgue measure (the density p = 1) and satisfy the conditions 
Z?2 and D' , sufficient to get extreme valuers distributions. For the second map, 
the algebraic automorphisms of the torus better known as the Arnold cat map, 
the existence of extreme value laws follows from the theory developed in [Gupta 



et al. 12009 . Starting from the definitions provided by Gnedenko we derive as 



a novel result the exact expression for the normalising sequences a m and b 7l 
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Case 1: gi(x)= -log(dist(x,£)). By equations [8] and [9] we know that: 



1 - F(u) = 1 - v(g(diat(x, ()) < u) 

= 1 - z/(- log(dist(x, C)) < u) 
= 1 - i/(dist(x,C) > e' u ) 



(13) 



and the last line is justified by using Lebesgue's Differentiation Theorem. 
Then, for maps with constant density measure, we can write: 



l-F{u)~v{B e -u{Q) = n d e 



— ud 



(14) 



where d is the dimension of the space and Qd is a constant. To use Gnedenko 
corollary it is necessary to calculate up 

up = sup{u; F(u) < 1} 

in this case up = +oo. 

Using Gnedenko equation [6] we can calculate G(t) as follows: 



1 - F(u) , [°° e~ ud , 1 [°° e~ v , 1 
" l " = ' ^ dlu ^ dV= d 



(15) 



According to the Leadbetter et al. 1983 proof of Gnedenko theorem we can 



study both a m and b m or j m convergence as: 

lim m(l - F{-f m + xG{j m )}) = e 

m— >oo 

lim mn d e- d ^ +xG ^ ) = e~ x 



(16) 



then we can use the connection between j m and normalising sequences to 
find a m and b m . 

By equation [5] or using relation |16| 



In 



so that: 

a m d 

Case 2: g2(x)=dist(x,C)~ 1 / a . We can proceed as for g\\ 



1 lnffirf) 
Om = j m(m) H — 



1 - F{u) = 1 - v(dist{x, CT 1,a < u) 
= 1 - i/(dist(af,C) > w" Q ) 

= KS„-<»(C)) = ^- ad 



(17) 
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in this case uf = +00. 



7m = F- 1 (l-l/m) = (mfi d ) 1 /M) (18) 

and, as discussed in section [2j using Beirlant |2004 choice of normalising 
sequences we expect: 



b„, = c ■ 



where c € K is a constant. 



Case 3: g3(x)=C-dist(x,C) 1 ^ Q; . Eventually we compute a m and b m for the g% 

observable class: 



1 - F(u) = l-v{C- dist(x, C) 1/Q < u) 
= l-i/(dist(x,C) > (C-u) a ) 
= v{B {c _ ur {0) = n d {C-uT d 

in this case uf = C . 



(19) 



7™ = - l/m) =C- (mfirf)- 1 ^ (20) 

For type 3 distribution: 

a m = {u F - 7™)~\ b m = u F ; (21) 
3.2 Numerical Experiments 

Since we want to show that unnormalised data may be fitted by using the GEV 
distribution Fq(x; ^t,cr, ^) we expect to find the following equivalence: 

a rn = l/cr b m = /! 

where, clearly, fi = fi(m) and a = o(m). This fact can be seen as a lin- 
ear change of variable: the variable y = a m (x — b m ) has a GEV distribution 
Fg(iJj A t = 0j cr = l>0 (that is an EV one parameter distribution with a m and b m 
normalising sequences) while x is GEV distributed Fq{x; [l = b m , a = l/a m , £). 

As we said above we now apply the previous considerations to two maps 
which enjoy extreme values laws and have constant density: we summarize 
below the theoretical results we obtained for all three type of observables. We 
have obtained the results in terms of m but, since we fix k = n ■ m, the previous 
results can be translated in terms of n as follows: 

For gi type observable: 

a = ^ A* °c 2 ln ( fc / n ) ( 22 ) 
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For (72 type observable: 



(23) 



For g 3 type observable: 



,l/(ad) 



fi = C (24) 

Following Freitas et al. 2009] we obtain the expression for the shape pa- 
rameters: £ = for gi type , £ = l/(ad) for g 2 type and £ = — l/(ad) for g 3 
type. 

In order to provide a numerical test of our results we consider a one-dimensional 
and a two dimensional map. The one dimensional map used is a Bernoulli Shift 
map: 



x t+i — 1 x t mod 1 5>leN 



(25) 



with q = 3. 



The considered two dimensional map is the famous Arnold's cat map defined 
on the 2-torus by: 







2 


f 




_yt+i_ 




i 


l 





mod 1 



(26) 



A wide description of properties of these maps can be found in Arnold and 



Avez 1968 and Hasselblatt and Katok 2003 



We proceed as follows. For each map we run a long simulation up to k it- 
erations starting from a given initial condition £. Note that the results - as we 
tested - do not depend on the choice of (. From the trajectory we compute the 
sequence of observables g\ , g 2 , g% as follows dividing it into n bins each contain- 
ing m — k/n observations. Then, we test the degree of agreement between the 
empirical distribution of the maxima and the GEV distribution according to the 
theoretical values presented above. A priori, it is reasonable to assume GEV 
as a suitable family of statistical models. For some selected values of n, the 
maxima are normalised and fitted to GEV distributions Fq(x; fi, rj, £) using a 
maximum likelihood method which selects values of the model parameters that 
produce the distribution most likely to have resulted in the observed data. 
All the numerical analysis contained in this work has been performed using 
MATLAB Statistics Toolbox functions such as gevfit and gevcdf. These func- 
tions return maximum likelihood estimates of the parameters for the generalized 
extreme value (GEV) distribution giving 95% confidence intervals for estimates 
Martinez and Martinez! 120021. 



As in every fitting procedure, it is necessary to test the a posteriori goodness 
of fit. We anticipate that in every case considered, fitted distributions passed, 
with maximum confidence interval, the Kolmogorov-Smirnov test described in 
For illustration purposes, we present in figure [I] an empirical 



Lilliefors 1967 
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pdf and cdf with the corresponding fits. 

Once k is set to a given value (in our case k — 10 7 ), the numerical simulations 
allow us to explore two limiting cases of great interest in applications where the 
statistical inference is intrinsically problematic: 

1. n is small (m is large), so that we extract only few maxima, each corre- 
sponding to a very extreme event. 

2. to is small (n is large), so that we extract many maxima but most of those 
will not be as extreme as in case 1). 

In case 1), we have only few data - of high quality - to fit our statistical 
models whereas in case 2) we have many data but the sampling may be spoiled 
by the inclusion of data not giving a good representation of extreme events. 
We have in general that in order to obtain a reliable fit for a distribution with 
p parameters we need 10 p independent data [Felici et al. 



2007a so that we 



expect that fit procedure gives reliable results for n > 10 . As the value of 
m determines to which extend the extracted bin maximum is representative of 
an extreme , below a certain value m m i n our selection procedure will be un- 
avoidably misleading. We have no obvious theoretical argument to define the 
value of m min . We expect to obtain good fits throughout the parametric region 
where the constraints on 77, m are satisfied. Therefore, our flexibility in choosing 
satisfying pairs (77, m) increases with larger values of k. 

For a g\ type observable function the behavior against 77 of the three param- 
eters is presented in figure [2] According to equation [22] we expect to find f = 0. 
For relatively small values of 77 the sample is too small to ensure a good con- 
vergence to analytical £ and confidence intervals are wide. On the other hand 
we see deviations from expected value as 777 < 10 3 that is when 77 > 10 4 . For 
the scale parameter a similar behavior is achieved and deviations from expected 
theoretical values a — 1/2 for Arnold Cat Map and a — 1 for Bernoulli Shift 
are found when 77 < 10 3 or m < 10 3 . Location parameter /i shows a logarithm 



decay with 77 as expected from equation 22 A linear fit of [i in respect to log(?7 



is shown with agray line in figure [2j The linear fit computed angular coefficients 



K* of equation 22 well approximate 1/d: for Bernoulli Shift map we obtain 
\K*\ = 1.001 ± 0.001 while for Arnold Cat map \K*\ = 0.489 ± 0.001. We find 
that £ values have best matching with theoretical ones with reliable confidence 
interval when both 77 > 10 3 and 777 > 10 3 . These results are confirmed even 
for g2 type and (73 type observable functions as shown in figures |3^) and [4^,) 
respectively. We present the fit results for a = 3 but we have done tests for 
different a and for fixed 77 and different a. 

For g 2 observable function we can also check that /1 and a parameters follow 



a power law as described in eq. 23 In the log-log plot in Figure |3p), |3p), we 
can see a very clear linear behavior. For the Bernoulli Shift map we obtain 
\K*\ = 0.330 ± 0.001 for n series , \K*\ = 0.341 ± 0.001 for a in good agree- 
ment with theoretical value of 1/3. For Arnold Cat map we expect to find 
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K* = 1/6, from the experimental data we obtain \K*\ = 0.163 ± 0.001 for \i 
and \K*\ =0.164 ±0.001 for a. 

Eventually, computing g^ as observable function we expect to find a constant 
value for /i while a has to grow with a power law in respect to n as expected 



from equation 24 As in g 2 case we expect \K*\ = 1 /(ad) and numerical 
results shown in figure |3|d) , |4j;) are consistent with the theoretical one since 
\K*\ = 0.323 ± 0.006 for Bernoulli shift map and \K*\ = 0.162 ± 0.006 for 
Arnold Cat map. 

In all cases considered the analytical behavior described in equation [23] and 
24 is achieved and the fit quality improves if n > 10 3 and m > 10 3 . The g% type 



observable constant has been chosen C = 10. The nature of these lower bound 
is quite different: 

4 Distributions of Extremes in mixing map with 
non-constant density measure 

4.1 Asymptotic sequences 

The main problem when dealing with maps that have absolutely continuous but 
non-constant density measure is in the computation of the integral: 

v(B s (0) = f P(x)dx (27) 

where Bs(() is the d-dimensional ball of radius S centered in (. 

We have to know the value of this integral in order to evaluate F(u) and, 

therefore, the sequences a m and b m . 

As shown in the previous section S is linked to the observable type: in all cases, 
since we substitute u = 1 — 1/m, 5 — > means that we are interested in m — > oo. 
In this limit, a first order approximation of the previous integral is: 

u(B s (0) * p(()5 d + 0(5 d+1 ) (28) 
that is valid if we are not in a neighborhood of a singular point of p(C)- 

As an example we compute the asymptotic sequences for a logistic map: 

x t+ i = rx t (l - x t ) (29) 
with r = 4. This map satisfies hypothesis described in the analysis per- 



formed for Benedicks-Carleson maps in Moreira Freitas and Freitas 1 2008 



For this map the density of the absolutely continuous invariant measure is 
explicit and reads: 

p(0 - J—; ? Ce(o,i) (30) 
7iVC(i - 
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So that: 



(31) 



2 r -I 

— arcsin(i/ £ + S — arcsin(-\/ £ — S 

where £ + 6 < 1 and £ — (5 > 0. Since Extreme Value Theory effectively 
works only if n,m are large enough, the results in eq. [31] can be replaced by a 
series expansion for S — > 0: 



2 



2d 



MVC + S ~ ^sm(^/C^S = [1 + J 2 P(C) + -] (32) 

7r vC(i-C) 



up to order <5 3 , where: 

p(c) = " c(w) + c 2 (i-o + c 2 (i-o 2 (33) 

Using the last two equations we are able to compute asymptotic normalising 
sequences a m and b m for all gi observables. 

Case 1: gi(x)= -log(dist(x,£)). For g\ observable functions we set S = e~ ud . 
In case of logistic map d — 1. First we have to compute G(t) using equation 15 
and the expansion in eq. : 



f°°du(e- u + e- 3u P(() 2 „ 

We can compute j m , if m » 1, as follows: 



F{ lm ) ~ 1 - - (35) 



At the first order in eq. [32] we get 

1 1 2e"T" 



m 71 VW- - 

so that: 



(36) 



1 ^ Hm) + ,I {WW = T) (37) 



Therefore, the sequences a m and 6 m if m >> 1 are: 



[G( 7m )]-^1+^ C (C- 1)P(0 (38) 



7 OT ~ln(m)+ln(2p(0) (39) 
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Case 2: g2(x)=dist(x,£) x l a . We can proceed as for g\ setting 5 — (au) a , 
computing j m we get at the first order in eq. |32| 

1 ( 1 Y 1/( " , ^ 

We can respectively compute a m and b m as: 

a m = 7™ b m = (2mp(C))-« (42) 

Case 3: g3(x)=C-dist(x,C) 1 / Q . As in the previous cases, we compute 7 m up 
to the first order setting S — [a{C — Jm)] a - 

If 1 \ 1/Q 

^ = c -«{^) (44) 

For type 3 distribution: 

a rn — ( U F — Tm)^ 1 ! b m = U F \ (45) 

where ut = C. 

4.2 Numerical experiment on the logistic map 

Following the same procedure detailed in section 3.2, we want to show the 
equivalence between EV computed normalising sequences a m and b m and the 
parameters of a GEV distribution obtained directly fitting the data even in case 
of logistic map that has not constant density measure. Using eq. [38p9] for g\ , 
we obtain the following theoretical expression: 

ff K()~l + ^((C-l)P(() /i(m,C)^ln(m)+ln(2p(C)) (46) 
From eq. [39] for g^ observable type, we write: 

<7(m,0^i(2mp(C))- MKC)-(2mp(C))^ (47) 
a 

and in g 3 case using eq. |45[ we expect to find: 

<T(m,0^-(2mp(C)r 1/a ti(m,()~C = u F (48) 
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Values of £ are independent on density and, as stated in Freitas' £ = for 
9i tyP e i £ = l/(otd) for g 2 type and £ = -l/(ad) for g 3 type. 

In figures [5][7] we presents a numerical test of the asymptotic behavior de- 
scribed in equations [46] - [48] on logistic map for d = 1 , a = 3, C = Up = 10, 
£ = 0.3 against the variable n. As shown in previous section, block maxima 
approach works well with maps with constant density measure when n and to 
are at least 10 3 : In fact, regarding £ parameter. Significant deviations from the 
theoretical value are achieved when n < 1000 or to < 1000 even in the case of 
the Logistic Map. 

Regarding /i and a, for 171 observable a linear fit of \i in respect to log(n) give 
us I If* I = 0.999 ± 0.002, while a shows the same behavior of £ since the best 
agreement with theoretical value a — 1 is achieved when n,m > 10 3 . In the log- 
log plots of figure ^jp) , |6p) for g 2 observable, we can observe again the expected 
linear behavior for /1 and a with \K*\ corresponding to l/(ad). From numerical 
fit we obtain \K*\ = 0.3334 ± 0.0007 for fj, series and \K*\ = 0.337 ± 0.002 for 
a in good agreement with theoretical value of 1/3. By applying a linear fit to 
the log-log plot in figure [Tfj) , the angular coefficient corresponding to a series 
is \K\ — 0.323 ± 0.003 again consistent with the theory. 

For a logistic map we can also check the GEV behavior in respect to initial 
conditions. If we fix n* = m* = 10 3 and fit our data to GEV distribution for 
10 3 different £ <E (0, 1) an asymptotic behavior is reached as shown from the 
previous analysis. For g\ observable function we have observed that the first 
order approximation works well for all three parameters. Deviation from this 
behavior are achieved for £ — ¥ 1 and £ — > as the measure become singular 
when we move to these points and we should take in account other terms of 
the series expansion. Numerically, we found that deviations from first order 
approximation are meaningful only if £ < 10~ 3 and £ > 1 — 10 -3 . Averaging 
over £ both £ and a we obtain < £ >= 1.000 ± 0.009 and < a >= 1.00 ± 0.03 
where the uncertainties are computed with respect to the estimator. Since 
we expect £ = and a = 1 at zero order approximation, numerical results 
are consistent with the theoretical ones; furthermore, experimental data are 
normally distributed around theoretical values. 

Asymptotic expansion also works well for g 2 observables: we obtain < £ >= 
0.334 ±0.001 in excellent agreement with theoretical value £ = 1/3. Eventually, 
in g 3 , averaging £ over different initial conditions we get < £ >= —0.334 ±0.002 
that is again consistent to theoretical value -1/3. 

5 The case of regular maps 

Freitas and Freitas 1 2008a] have posed the problem of dependent extreme values 
in dynamical systems that show uniform quasi periodic motion. Here we try 
to investigate this problem numerically. We have used a one-dimensional and a 
bi-dimensional discrete map. The first one is the irrational translation on the 
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torus defined by: 



x t+ i =x t + /3 mod 1 /3 e [0, 1] \ Q 
And for the tridimensional case, we use the so called standard map: 



(49) 



Vt+i = Vt + 7T sin(27rx t ) mod 1; x t +i = x t + y t+ i mod 1. (50) 
Air 

with A = 10~ 4 . For this value of A, the standard map exhibits a regular 
behavior and it is not mixing, as well as torus translations. This means that 
these maps fail in satisfying hypothesis D2 and D' and moreover they do not 
enjoy as well an exponential hitting time statistics. About this latter statistics, 
it is however known that it exists for torus translation and it is given by a par- 
ticular piecewise linear function or a uniform distribution depending on which 



sequence of sets Ay. is considered Coelho and De Faria 1996 . In a similar way, 



a non-exponential Hitting Time Statistics (HTS) is achieved for standard map 
when A << 1 as well as for a skew map, that is a standard map with A = 
|Buric et al. 2005 . Therefore we expect not to obtain a GEV distribution of 
any type using gi observables. 

We have pointed out that the observable functions choice is crucial in order 
to observe some kind of distribution of extreme values when we are dealing 



with dynamical systems instead of stochastic series. Nicolis et al. 2006 have 



shown how it is possible to obtain an analytical EV distribution which does not 
belong to GEV family choosing a simple observable: they considered the series 
of distances between the iterated trajectory and the initial condition. Using the 
same notation of section [2] we can write: 

Y m (x = fC) = dist(/*C, C) Mm = mm{Y , -Y m -i} 

For this observable they have shown that the cumulative distribution F(x) — 
P{a m (M m — b m ) < x} of a uniform quasi periodic motions is not smooth but 
piecewise linear (Nicolis et al. 1 2006 , figure 3). Furthermore slop changes of 
F(x) can be explained by constructing the intersections between different iter- 
ates of equation 49 F(x) must correspond to a density distribution continuous 
obtained as a composition of box functions: each box must be related to a 
change in the slope of F(x). 

The numerical results we report below confirm that for the maps [49] and 
1501 the distributions of maxima for various observables cannot be fitted with a 
GEV since they are multi modal. We recall that the return times into a sphere 
of vanishing radius do not have a spectrum, if the orbits have the same fre- 
quency, whereas a spectrum appears if the frequency varies continuously with 



the action, as in the standard map for A close to zero |Hu et al. 2004 . Since 
the EV statistics refers to a single orbit, no change due to the local mixing, 
which insures the existence of a return times spectrum |Hu et ah 2004 , can be 
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observed. Considering that the GEV exists when the system is mixing and does 
not when it is integrable, one might use the quality of fit to GEV as a dynamical 
indicator, for systems which exhibit regions with different dynamical properties, 
ranging from integrable to mixing as it occurs for the standard map when A is 
order 1. Indeed we expect that in the neighborhood of a low order resonance, 
where the omoclinic tangle of intersecting separatrices appears, a GEV fit is 
possible. Preliminary computations carried out for the standard map and for 
a model with parametric resonance confirm this claim, that will be carefully 
tested in the near future. 



Using the theoretical framework provided in Nicolis et al. |2006 we check 
numerically the behavior of maps described in eq. 49j|50 analysing EV distribu- 
tions for gi observable functions. Proceeding as in section 3 for mixing maps, 
we try to perform a fit to GEV distribution starting with different initial condi- 
tions C) a set of different a values and (n, m) combinations. In all cases analysed 
the Kolmogorov Smirnov test fails and this means that GEV distribution is not 
useful to describe the behavior of this kind of statistics. This result is in agree- 
ment with Freitas et al. 2010b but we may find out which kind of empirical 
distribution is obtained. 



Looking in details at M m histograms that correspond to empirical density 
distributions, they appear always to be multi modal and each mode have a well 
defined shape: for g\ type observable function modes are exponential while, for 
g2 and 173, their shape depends on a value of observable function. Furthermore, 
the number of modes and their positions are highly dependent on both n, in and 
initial conditions. 

Using Nicolis et al. 2006 results it is possible to understand why we obtain this 
kind of histograms: since density distribution of M m is a composition of box 
functions, when we apply observables we modulate it changing the shape of 
the boxes. Therefore, we obtain a multi modal distribution modified according 
to the observable functions gi. 

An example is shown in figure [8] for standard map: the left figures correspond 
to the histogram of the minimum distance obtained without computing gi ob- 
servable and reproduce a composition of box functions. The figures in the right 
show how this distribution is modified by applying g\ observable to the series 
of minimum distances. We can see two exponential modes, while the third is 
hidden in the linear scale but can be highlighted using a log-scale. The up- 
per figures are drawn using n = 3300, m — 3300, the lower with n = 10000, 
m = 1000. 



6 Concluding Remarks 

EVT was developed to study a wide class of problems of great interest in different 
disciplines: the need of modeling events that occur with very small probability 
comes from the fact that they can affect in a strong way several socioeconomic 
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activities: floods, insurance losses, earthquakes, catastrophes. A very extensive 
account of EVT applications has been recently given in Ghil et al. [2011 . EVT 
was applied on limited data series using the block-maxima approach facing the 
problem of having a good statistics of extreme values retaining a sufficient num- 
ber of observation in each bin. Often, since no theoretical a priori values of 
GEV parameters are available for this kind of applications, we may obtain a 
biased fit to GEV distribution even if tests of statistical significance succeed. 
The recent development of an extreme value theory in dynamical systems give 
us the theoretical framework to test the consistency of block-maxima approach 
when analytical results for distribution parameters are available. This theory 
relies on the global properties of the dynamical systems considered (such as the 
degree of mixing or the decay rate of the Hitting Time Statistics) but also on 
the observable functions we chose. 



Our main finding is that a block-maxima approach for GEV distribution 
is totally equivalent to fit an EV distribution after normalising sequences are 
computed. To prove this we have derived analytical expressions for a rn and b rn 
normalising sequences, showing that fi and a of fitted GEV distribution can 
replace them. This approach works for maps that have an absolutely continu- 
ous invariant measure and retain some mixing properties that can be directly 
related to the exponential decay of HTS. Since GEV approach does not require 
the a-priori knowledge of the measure density that is instead require by the EV 
approach, it is possible to use it in many numerical applications. 

Furthermore, if we compare analytical and numerical results we can study 
what is the minimum number of maxima and how big the set of observations 
in which the maximum is taken has to be. To accomplish this goal we have 
analysed maps with constant density measure finding that a good agreement 
between numerical and analytical value is achieved when both the number of 
maxima n and the observations per bin m are at least 10 3 . We remark that 
the fits have passed Kolmogorov Smirnov test with maximum confidence interval 
even if n < 10 3 or < m < 10 3 so that parametric or non parametric tests are not 
the only thing to take in account when dealing with extreme value distributions: 
if maxima are not proper extreme values (which means m is not large enough) 
the fit is good but parameters are different from expected values. The lower 
bound of n can be explained using the argument that a fit to a 3-parameters 
distribution needs at least 10 3 independent data to give reliable informations. 
Therefore, we checked that in case of non-constant absolutely continuous den- 
sity measure the asymptotic expressions used to compute \i and a works when 
we consider n and m of order 10 3 . For logistic map the numerical values of 
parameters we obtain averaging over different initial conditions are totally in 
agreement with the theoretical ones. In regular maps, as expected, the fit to a 
GEV distribution is unreliable. We obtain a multi modal distribution, that, for 
the analyzed maps, is the result of a composition of modes in which the shape 
depends on observable types. This behavior can be explained pointing out that 
this kind of systems have not an exponential HTS decay and therefore have no 
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EV law for observables of type gi . 

To conclude, we claim that we have provided a reliable way to investigate 
properties of extreme values in mixing dynamical systems which may satisfy 
mixing conditions (like D 2 and D'), finding an equivalence among a m , b m , [i 
and o~ behavior for absolutely continuous measures. In our future work we intend 
to address the case of singular measure. Recently the theorem was generalised 
to the case of non smooth observations and therefore it holds also with non ab- 



solutely continuous invariant probability measure Ereitas et al. 2010b . In this 
case we expect the same for all the procedure described here. Understanding 
the extreme values behavior for singular measures will be crucial to apply pro- 
ficiently this analysis to operative geophysical models since in these case we are 
always dealing with singular measures. In this way we will provide a complete 
tool to study extreme events in complex dynamical systems used in geophysical 
or financial applications. 
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Fig. 1: Left: g\ observable empirical histogram and fitted GEV pdf. Right: g\ 
observable empirical cdf and fitted GEV cdf. Logistic map, n = 10 4 , 
to = 10 4 
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Fig. 2: <?i observable, C ^ 0.51. a) f VS log 10 (n); b) cr VS log 10 (n); c) /x VS 
log(n). Right: Bernoulli Shift map. Left: Arnold Cat Map. Dotted lines 
represent computed confidence interval, gray lines represent linear fits 
and theoretical values. 
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Fig. 3: 32 observable, C - 0.51. a) £ VS log 10 (n); b) log 10 (<r) VS log 10 (n); 
c) log 10 (/i) VS log 10 (n). Right: Bernoulli Shift map. Left: Arnold Cat 
Map. Dotted lines represent computed confidence interval, gray lines 
represent linear fits and theoretical values. 



7 Acknowledgments 



23 




Fig. 4: 33 observable, C - 0.51. a) £ VS log 10 (n); b) log 10 (<r) VS log 10 (n); 
c) log 10 (/i) VS log 10 (n). Right: Bernoulli Shift map. Left: Arnold Cat 
Map. Dotted lines represent computed confidence interval, gray lines 
represent linear fits and theoretical values. 



7 Acknowledgments 



24 




Fig. 5: g\ observable, £ = 0.31. a) £ VS log 10 (n); b) a VS log 10 (n); c) /i 
VS log(n). Logistic map. Dotted lines represent computed confidence 
interval, gray lines represent linear fits and theoretical values. 
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Fig. 6: g 2 observable, ( = 0.3. a) £ VS log 10 (ra); b) log 10 (cr) VS log 10 (n); c) 
log 10 (/i) VS log 10 (n). Logistic map. Dotted lines represent computed 
confidence interval, gray lines represent linear fits and theoretical values. 
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Fig. 7: g 3 observable, ( = 0.3. a) £ VS log 10 (ra); b) log 10 (cr) VS log 10 (n); c) 
log 10 (/i) VS log 10 (n). Logistic map. Dotted lines represent computed 
confidence interval, gray lines represent linear fits and theoretical values. 
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a) 




0.2 0.4 0.6 0.8 1 1.2 5 10 15 20 

min(dist(x,Q) x 1Q -3 ^= -l,og(min(dist(x,Q)) 



b) 




0.5 1 1.5 5 10 15 20 

min(dist(x,Q) x 1Q -3 g^ -log(min(dist(x,Q)) 



Fig. 8: Histogram of maxima for g\ type observable function, standard map, 
x o = Ho = v2 — 1- Left: series of min(dist(/ t C, C))- Right: series of 
,9i = -log(min(dist(/*C,C)))- a) n = 3300, m = 3300. b) n = 10000, 
m = 1000. 
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